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ABSTRACT 

Discerning the likelihood of the so-called runaway instability of thick accre¬ 
tion disks orbiting black holes is an important issue for most models of cosmic 
gamma-ray bursts. To this aim we investigate this phenomenon by means of 
time-dependent, hydrodynamical simulations of black hole plus torus systems in 
general relativity. The evolution of the central black hole is assumed to be that 
of a sequence of Kerr black holes of increasing mass and spin, whose growth rate 
is controlled by the transfer of mass and angular momentum from the material 
of the disk spiralling in through the event horizon of the black hole. The self¬ 
gravity of the disk is neglected. We find that when the black hole mass and spin 
are allowed to increase, constant angular momentum disks undergo a runaway 
instability on a dynamical timescale (a few orbital periods). However, our simu¬ 
lations show that a slight increase of the specihc angular momentum of the disk 
outwards has a dramatic stabilizing effect. Our results, obtained in the frame¬ 
work of general relativity, are in broad agreement with earlier studies based both 
on stationary models and on time-dependent simulations with Newtonian and 
pseudo-Newtonian gravitational potentials. 

Subject headings: Accretion — Accretion disks — Black hole physics — Gamma- 
rays: bursts — Hydrodynamics — Relativity 


1. Introduction 


The so-called runaway instability of thick accretion disks orbiting black holes was first 
noticed by Abramowicz et ah (1983). The origin of the instability is a dynamical process 
by which the cusp of a disk initially hlling its Roche lobe moves outwards due to mass 
transfer from the disk to the accreting black hole. This process leads to the complete 
destruction of the disk on a dynamical timescale. Abramowicz et al. (1983) found that 
the runaway instability occurs for a large range of parameters such as the disk-to-hole mass 
ratio, Md/Mbh, and the location of the disk inner radius. More detailed studies followed, 
most of which are based on stationary models in which a fraction of the mass and angular 
momentum of the initial disk is transferred to the black hole. The new gravitational held 
allows to compute the new position of the cusp, which controls the occurrence of the runaway 
instability. 

The main conclusion of those studies is that the very existence of the instability re¬ 
mains uncertain. Taking into account the self-gravity of the disk seems to favor the insta¬ 
bility, as shown from both, studies based on a pseudo-potential for the black hole (Khanna 
& Chakrabarti 1992; Masuda et al. 1998) and from fully relativistic calculations (Nishida 
et al. 1996). However, the rotation of the black hole has a stabilizing effect (Wilson 
1984; Abramowicz et al. 1998), and the same happens when the disks are built with non¬ 
constant distributions of the specihc angular momentum (increasing outwards) (Daigne & 
Mochkovitch 1997; Abramowicz et al. 1998). 

In a recent paper (Font & Daigne 2002) we presented the hrst time-dependent, hydro- 
dynamical simulations in general relativity of the runaway instability of constant angular 
momentum thick disks around a Schwarzschild black hole. The self-gravity of the disk was 
neglected and the evolution of the central black hole was assumed to follow a sequence of 
Schwarzschild black holes of increasing mass. We found that by allowing the mass of the 
black hole to grow the runaway instability appears on a dynamical timescale, in agreement 
with previous estimates from stationary models. For a black hole of 2.5 Mq and disk-to-hole 
mass ratios between 1 and 0.05, our simulations showed that the timescale of the instability 
never exceeds 1 s for a large range of mass fluxes and it is typically a few 10 ms. Similar 
results for different initial data have been recently reported by Zanotti et al. (2002). 

In this Letter we extend our study to the most interesting case of non-constant angular 
momentum disks. The main motivation of our work is to check through time-dependent 
simulations in general relativity whether such distributions have indeed the stabilizing effect 
previously found in non self-gravitating stationary models (Daigne & Mochkovitch 1997; 
Abramowicz et al. 1998; Lu et al. 2000). 
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2. Numerical framework 

A given initial state of a thick disk orbiting a black hole is determined by hve parameters: 
the black hole mass Mbh, the specihc angular momentum in the equatorial plane /, the 
potential barrier at the inner edge of the disk Ahhin = hhin — hhcusp, the polytropic constant n 
and the adiabatic index 7 of the equation of state (EoS). We assume a specihc law of rotation 
of the disk in which the angular momentum in the equatorial plane increases outwards as a 
positive power law / oc r". The specihc values we have considered for the exponent a are 
listed in Table 1. The initial equilibrium conhgurations are built in various steps (see Daigne 
& Font ( 2002 ) for details): hrst, from the coefficients of the Schwarzschild metric and the 
radial prohle of I in the equatorial plane, the angular momentum distribution outside the 
equatorial plane is computed by solving the equations describing the constant I surfaces, i.e. 
the von Zeipel cylinders. Then, the total specihc energy ut is obtained from the distribution 
of / via the normalization of the 4-velocity, u^j,u^ = — 1. The integral form of the relativistic 
Euler equation allows then to compute the quantity ^ = — dP/w, where P is the pressure 
and w the enthalpy. The Newtonian limit of W is the total (gravitational plus centrifugal) 
potential. The density and the pressure can be easily computed from W using the EoS. The 
geometrical structure of the equipotential surfaces is comparable to that of the Roche lobes 
of a binary system. In particular, there is a cusp where mass transfer from the disk to the 
black hole is possible. 

As we did in our previous investigation (Font & Daigne 2002) we focus on models which 
are expected to exist at the central engine of gamma-ray bursts (GRBs), formed either after 
the coalescence of a compact binary system (composed of two neutron stars or a neutron 
star and a stellar-mass black hole), or after the gravitational collapse of a massive star (see 
e.g. Woosley (2001)). Taking into account the various results from numerical simulations of 
those systems (Ruffert & Janka 1999; Kluzniak & Lee 1998; MacFadyen & Woosley 1999; 
Shibata & IJryu 2000 ; Aloy et al. 2000 ) we hx the mass of the black hole to Mbh = 2.5 M©. 
The angular momentum I is adjusted to yield a disk-to-hole mass ratio of 1 (models A and 


Table 1: Model parameters 


Model 

Vbh 

IM0I 

Md/Mbh 

a 


ffistat 

|Me/s] 

'^cusp 

^ center 

^orb 

[ms] 

trnn/torh 

A 

2.5 

1.0 

0 . 

3.81 X 10-2 

~ 25. 

4.896 

7.592 

1.62 


B 

2.5 

1.0 

0.1 

2.26 X 10-2 

~ 2.5 

5.225 

9.982 

2.44 

stable 

C 

2.5 

0.1 

0 . 

2.60 X 10-2 

~ 5.9 

4.962 

7.459 

1.58 

~ 8 

D 

2.5 

0.1 

0.075 

1.87 X 10-2 

~ 1.3 

5.243 

8.964 

2.08 

stable 



B) or of 0.1 (models C and D). The latter case is more realistic, being the parameters closer 
to those inferred from simulations of compact binary mergers. Furthermore, neglecting the 
disk self-gravity is also much better justihed for Md/Mbh = 0.1. An even more realistic 
model would assume an initially rotating black hole, as both in collapsars and mergers the 
hole is expected to form with an initial spin parameter > 0.5. Such models will be considered 
in a forthcoming study. 

The EoS adiabatic index and polytropic constant are, respectively, 7 = 4/3 and k, = 
4.76 X 10^^ cgs. This corresponds to an EoS dominated by the contribution of relativistic 
degenerate electrons (the typical density in the disk is ~ lO^^-lO^^ g/cm^). The gravitational 
potential at the inner edge of the disk is hxed to AfFin = 0.75|lTcusp| in models A and B, 
and 0.5|lFcusp| in models C and D, which corresponds to a hnite disk (IFin < 0), initially 
overflowing its Roche lobe (AfFin > 0) so that mass transfer is possible at the cusp. The 
parameters describing the models are summarized in Table 1. 

The equilibrium initial data are evolved in time using the same relativistic hydrodynam¬ 
ics code employed in Font & Daigne (2002). This code integrates the non-linear hyperbolic 
system of the general relativistic hydrodynamic equations in a Kerr spacetime, using Boyer- 
Lindquist (t,r,9,(j)) coordinates and under the restriction of axisymmetry, = 0. The nu¬ 
merical scheme, based on an approximate Riemann solver, is second-order accurate in both 
space and time due to the use of a piecewise linear reconstruction algorithm and a conser¬ 
vative, two-step Runge-Kutta scheme for the time update. The simulations use a numerical 
grid of 400 x 100 zones in r and 6, respectively. The radial grid is logarithmically spaced to 
account for sufficient resolution near the horizon (rmin = 2.12 Mbh, Armin = 1.99 x 10 “^Mbh; 
G = c = 1) and includes the whole disk. 

The procedure we follow to account for the spacetime dynamics is equivalent to the 
one employed in Font & Daigne (2002), with the important difference that now, both the 
black hole mass and its angular momentum are allowed to increase. This leads to additional 
metric coefficients - the 0 component of the shift vector - to be non zero, making the whole 
numerical treatment appreciably more difficult than in the case of the Schwarzschild metric. 
The spacetime metric is hence approximated at each time step by a stationary (Kerr) black 
hole metric of varying mass Mbh and angular momentum Jbh- We note that both quantities 
increase very slowly during the evolution. Details of our procedure as well as code tests for 
rotating black holes will be reported in Daigne & Font (2002). 



- 5 - 

3. Results 

Each of our four initial models is evolved three times. In the hrst series of runs Mbh 
and Jbh are kept constant, hence the spacetime is held hxed. In the second series of runs 
only Mbh is allowed to increase while in the third series of runs both, Mbh and Jbh are 
allowed to increase. The growth rate is monitored through the mass and angular momentum 
transfer from the disk to the black hole across the innermost grid zone. Since in realistic 
disks angular momentum is transported outwards by dissipative processes or removed from 
the system by gravitational radiation, we adopt a conservative value for the efficiency of 
the angular momentum transfer, assuming that only 20% of the angular momentum of the 
accreted material is transferred to the black hole. Results for other values of the efficiency 
are reported in Daigne & Font (2002). 

The evolution of the mass flux m = 27i y^—gDv^d9 for all twelve simulations is plotted 
in Fig. 1. In this expression g is the determinant of the Kerr metric, is the radial component 
of the 3-velocity and D = pT is the relativistic mass density, p being the rest-mass density 
and T the Lorentz factor. At early times the mass flux evolution for all three series in each 
model is exactly the same, irrespective of the increase of the black hole mass and spin being 
taken into account or not. For the hrst series of runs the mass hux reaches a stationary 
regime where m oc AhFd (Font & Daigne 2002), which is the theoretical expectation for 
7 = 4/3 (Kozlowski et al. 1978). Since we consider only models with hxed ratios M^/Mbh 
and AlTin/lTcusp, the stationary mass hux is diherent for each model. However our results 
are not modihed when all models have the same initial m (see Daigne & Font (2002)). 

For constant angular momentum disks (models A and C; a = 0) the time evolution of 
the system changes dramatically when the spacetime dynamics is taken into account, i.e. 
when either Mbh only or Mbh Jbh increase. In either case the runaway instability, 
rehected in the rapid growth of the mass accretion rate, appears on a dynamical timescale 
of ~ 4 forb ~ 6.5 ms for model A (~ 8 forb ~ 12 ms for model C). The time derivative of the 
mass hux also increases, which implies the rapid divergence of rh. Comparing the case where 
only Mbh changes with the case where both Mbh and Jbh change, the instability appears 
slightly later in the latter case (dotted line in Fig. 1), the qualitative behaviour of the mass 
hux being, however, identical. 

The main diherences appear in the case of non-constant angular momentum disks (mod¬ 
els B and D). In Fig. 1 it becomes apparent that a slight increase of the specihc angular 
momentum outwards (well below the Keplerian limit a = 0.5) strongly stabilizes the disk and 
completely suppresses the runaway instability. During the accretion process, the material 
at the inner edge of the disk is replaced by new material with a higher angular momentum. 
Therefore, the centrifugal barrier which acts against accretion is much more efficient. We 
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note that the long-term behaviour of the mass flux evolution simply reflects the fact that 
the initial accretion rate is non zero (highly super-Eddington: ~ 2.5 Mq/s for model B, 
~ 1.3 Mq/s for model D). Hence, by integrating m along the entire time of the simulation, 
it is possible to check that the total mass transferred becomes, asymptotically, of the order 
of the initial mass of the disk, which results in the end of the accretion process. This effect 
is more pronounced for smaller disk-to-hole mass ratios (model D). Notice that the mass 
flux decreases during the long-term evolution. This is due to the fact that for non-constant 
angular momentum disks (models B and D), the accretion process makes the cusp move 
faster towards the black hole than the inner radius of the disk (which is precisely why the 
instability is physically suppressed). Therefore, the potential at the inner edge AfTin de¬ 
creases with time and, consequently, also the mass flux. For this reason, accretion becomes 
increasingly slower, which makes numerically difficult to follow the evolution until the entire 
disk has been fully accreted. 

In Fig. 2 we plot the time evolution of the mass of the disk for all twelve simulations. 
The unstable behaviour of models A and C, with a = 0, is characterized by the sudden 
decrease of in less than ~ 6 — 10 orbital periods. By the end of the simulation the disk 
and black hole masses are ~ 0.13 Mq and Mbh ~ 4.88 Mq for model A (M^ ~ 0.013 
Mq and Mbh ~ 2.74 Mq for model C), which means that 95% of the disk mass has already 
been accreted. However, when the spacetime is held hxed (solid line), the mass of the disk 
for model A at t ~ 200 torb is still Md ~ 0.1 Mq (correspondingly, M^ ~ 0.015 Mq at 
t ~ 100 torb for niodel C). This means that the accretion timescale is at least 30 and 10 
times longer, for models A and C respectively, than the runaway (dynamical) timescale trun- 
On the other hand, for models B (a = 0.1) and D (a = 0.075), the stability properties of the 
non-constant distribution of angular momentum are dramatically reflected in the small loss 
of the mass of the disk throughout the evolution. At the end of the simulation the mass loss 
is only ~ 2% for model B and ~ 14% for model D (the lines are almost indistinguishable). 
The hnal masses of the black holes are Mbh ~ 2.55 Mq and ~ 2.54 Mq for models B and 
D, respectively. 

The transfer of mass and angular momentum from the disk to the black hole leads to 
the gradual increase of the rotation law of the initially non-rotating black hole. In Fig. 3 
we plot the time evolution of the black hole spin a = Jbh/M|jj for both the stable and 
unstable cases. The stable disks, models B and D, only transfer a very small fraction of their 
angular momentum to the black hole. At the end of the simulation the black hole spin is 
a = 0.014 for model B and a ~ 0.005 for model D. On the contrary, for the unstable disks 
the transfer of angular momentum is considerably higher, accelerating rapidly during the 
runaway instability. Therefore, in ~ 6 — 10 orbital periods the initial Schwarzschild black 
hole turns into a mildly (slowly) rotating Kerr black hole with a ~ 0.4 (a ~ 0.05) for model 
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A (model C). 


4. Discussion 

By means of time-dependent hydrodynamical simnlations in general relativity we have 
shown that the rnnaway instability of (non self-gravitating) thick disks around (rotating) 
black holes (Abramowicz et al. 1983) is strongly suppressed when the distribution of angular 
momentum in the disk is non constant. A small increase outwards of the specihc angular 
momentum distribution, well below the Keplerian limit, could result in a strong stabilizing 
effect. In particular, we have demonstrated that for disks whose angular momentum follows 
a power-law, I oc r", the case a = 0 is unstable for both, Md/Mbh = 1 and 0.1, while 
disks with small non-zero values of a are stable {a = 0.1 for Md/Mbh = 1 and a = 
0.075 for Md/Mbh = 0.1). This result is in complete agreement with previous studies 
based on stationary sequences of equilibrium conhgurations, either using a pseudo-Newtonian 
potential (Daigne & Mochkovitch 1997) or in general relativity (Abramowicz et al. 1998). 
For the hrst time we have shown this effect with a time-dependent relativistic calculation, 
and we have been able to estimate the lifetime of a thick disk orbiting a black hole in both, 
the unstable and stable cases. 

According to our results the runaway instability is most likely avoided in non rigidly 
rotating disks formed in the coalescence of a binary neutron star system or in the gravitational 
collapse of a massive star. If the growth of non-axisymmetric modes in the disk by the 
Papaloizou-Pringle instability is suppressed by the accretion process itself, as suggested by 
Blaes (1987), systems consisting of a Kerr black hole surrounded by a high density torus 
may then be long lived. The lifetime is probably controlled by the viscous timescale (a few 
seconds) rather than by the dynamical one, which may provide enough time for any plausible 
magneto-hydrodynamical process to efficiently transfer part of the energy reservoir of the 
system to a relativistic outflow. Therefore, the most favoured current GRB models (see 
Woosley (2001)) can indeed be based on such central engines. This is especially important 
in the context of the so-called internal shock model for the prompt gamma-ray emission (Rees 
& Meszaros 1994), where the observed lightcurve reflects the activity of the central engine, 
with no modihcation of the timescales other than the time dilation due to the redshift. In 
particular, the central engine has to survive for a duration at least comparable with the 
observed duration of the GRB. 

We note that our study does not include the self-gravity of the disk, despite its rele¬ 
vance when Md/Mbh ~ 1- Successful attempts to long-term stable simulations of thick disks 
orbiting black holes by solving the coupled system of Einstein and general relativistic hydro- 



dynamic equations seem, however, out of the scope of present day numerical relativity codes 
(see, e.g. Font (2000) and references there in). According to previous studies (Khanna & 
Chakrabarti 1992; Nishida et ah 1996; Masuda et ah 1998), based on either time-dependent 
simulations with pseudo-potentials or stationary models in general relativity, self-gravity has 
a strong destabilizing effect. However, this effect has to decrease for small disk-to-hole mass 
ratios. Therefore, we expect that for Md/Mbh ^ 0.1 the stabilizing effect of non-constant 
angular momentum distribution demonstrated in this study could overcome the destabilizing 
effect of the much less relevant disk’s self-gravity. 

A comprehensive study of a broader sample of initial models accounting for different 
disk-to-hole mass ratios and initial accretion rates will be presented elsewhere (Daigne & 
Font 2002). We further plan to perform a careful study of the critical value of a separating 
stable and unstable disks and to hnd its dependence on the disk-to-hole mass ratio and on 
the black hole spin. 
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Fig. 1.— Time evolution of the mass flux. The solid lines correspond to evolutions where 
the spacetime is held hxed, while the dashed and dotted lines correspond to the cases where 
only the growth of Mbh is allowed, and where both Mbh and Jbh increase, respectively. 
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Fig. 2.— Time evolution of the mass of the disk. The line style follows the same convention 
as in Fig. 1. The non-constant angnlar momentnm disks, models B and D, remain stable 
throughont the whole simulation. 
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Fig. 3.— Time evolution of the spin of the black hole for the simulations in which the 
transfer of mass and angular momentum from the disk is allowed. 



